Load data

Data is saved as an excel file. This is a fine example of dirty messy data created during research.

messy

This contains data from multiple experiments (causing batch effects), it has missing data for some variables. It has variables whose name does not convey much meaning. This includes the variables ‘OK for strain comparison’ and ‘remove from analysis (justified)’. These are QC columns that would not appear in publication quality datasets but often appear in raw datasets. In this respect we are adhering to a core principle of data science never touch the raw data. In truth, this is not really raw data, it is a compilation of raw data (cell counts, worm counts, mouse IDs) and analysed data (produced in flow cytometry analysis) painstakingly complied into a curated dataset. However because it was is precious in it’s current form we will treat the excel file as sacrosanct - do not edit.

To make things simple for this course, I have deleted most of flow cytometry data leaving only the % proportion of each cell type.

Now lets load the data. This will build on what you learned in course 1. There are too many columns (cols) to parse each one individually. The read_excel() function will automatically parse the data and assign it to the correct data type (we hope). Much easier than specifying each col. 

#make sure you don't have the file open in excel!
raw <- read_excel("./data/plueral_cavity_database.xlsx")

What size is the tibble? Lets use the dim() or dimension function. Btw, a tibble is the tidyverse version of a data frame.

dim(raw) #rows x cols

[1] 666  29

Each row represents a mouse and each col an individual datapoint about that mouse. Lets print the colnames with the automatically assigned datatype by the read_excel() function.

print(sapply(raw,class)) # this is an example of a vectorised function use in R 

                       Global ID                       Experiment 
                       "numeric"                      "character" 
                           Mouse               Experimental Group 
                       "numeric"                      "character" 
                       Infection                              Sex 
                     "character"                      "character" 
                          Strain                       Background 
                     "character"                      "character" 
                        Genotype                        Treatment 
                     "character"                      "character" 
        OK for strain comparison                        Timepoint 
                     "character"                        "numeric" 
Remove from Analysis (justified)                        Day 0 Age 
                     "character"                        "numeric" 
                             Age                       Live worms 
                       "numeric"                        "numeric" 
                      Cell Count                          B cells 
                       "numeric"                        "numeric" 
                     Eosinophils                      Neutrophils 
                       "numeric"                        "numeric" 
                    CD4+ T cells                     CD4- T cells 
                       "numeric"                        "numeric" 
                         T cells                               DC 
                       "numeric"                        "numeric" 
                   Myeloid cells                        Monocytes 
                       "numeric"                        "numeric" 
                             SCM                   Converting SCM 
                       "numeric"                        "numeric" 
                             LCM 
                       "numeric" 

Its too long to print inline but you can look at it in a spreadsheet format using view(raw). Looking at the col names we can see that some of the cols are character vectors (categorical). These include cols that are likely to be metadata, such as ‘Experiment’ (which experiment the datapoint comes from). Other character variables could be considered to be results, ‘sex’ and ‘age’. Most of the variable are numeric. The numeric data includes cell counts, worm counts and flow cytometry data. The flow cytometry data was generated using the following gating, and expressed as % of total live cells.

Create categories

Factorising

Factors in R are basically a way of assigning categories to the data, so instead of being stored as a new character string for each row, they are given a unique ID for each factor that is used repeatedly. This has lots of benefits, not least it is computationally cheaper requiring less memory.

# see what the options are.
table(raw$Infection)

Infected    Naive 
     313      161 

table(raw$Background)

 BALB/c C57BL/6 
    201     273 

if we look at the current class of the variable it is a character vector and has no factor ‘levels’

class(raw$Infection)

[1] "character"

levels(raw$Infection)

NULL

Now we will factorize and check the class of the updated variable

raw$Infection <- factor(raw$Infection,levels = c("Naive", "Infected")) # notice we define the order 
raw$Background <- factor(raw$Background,levels = c("C57BL/6", "BALB/c"))
class(raw$Infection)

[1] "factor"

We can investigate the facrors using the levels() function

levels(raw$Infection) # we could also run this like this : raw$Infection %>% levels()

[1] "Naive"    "Infected"

Basic plotting

We will do this by adding layers to the ggplot using the + operator. Let’s start with the basic ggplot function. Here we pass the dataset to ggpot but not told ggplot to do anything, so all we get is a blank plot area.

raw %>% ggplot(aes())

Next we tell ggplot to plot infection status (as a categorical variable) on the x axis and % myeloid cells on the y axis. For this we use the aes() within the ggplot function. However while ggplot shows those axies it doesn’t show any data. That is because ggplot has no default way it displays data, we must tell ggplot to plot a bar chart or plot points, etc.

raw %>% ggplot(aes(x = Infection, y=`Myeloid cells`))

Adding points

Below we tell ggplot to plot points using the geom_point() function. Note that we use an + to add additional elements to the graph, these are called layers. Finally we see the data. It’s not very attractive however, with the points all bunched up and we can’t really infer anything about data.

raw %>% ggplot(aes(x = Infection, y=`Myeloid cells`)) +  
  geom_point()

boxplots

Boxplots are an excellent way to view data which have high n number. We add these using geom_boxplot()

raw %>% ggplot(aes(x = Infection, y=`Myeloid cells`)) + 
  geom_point() +
  geom_boxplot()

Next we want to add colour to the plot. To do this globally (i.e. to all layers) we add it within the aes() function as aes(col=Background). This colours both our boxplots, points and legend by genetic background.

raw %>% ggplot(aes(x = Infection, y=`Myeloid cells`, col=Background)) + 
  geom_point() +
  geom_boxplot()

fancy points

The col argument actually separated out boxplot but not the points. For this we will make use of a new R package called ggbeeswarm. You can do this without an additional package but I like this method as the way the points get distributed is very similar to the way Prism handles points.

raw %>% ggplot(aes(x=Infection, y= `Myeloid cells`, col=Background)) +
  geom_boxplot() + 
  geom_quasirandom(method = "tukeyDense", dodge.width = 0.75)

Much nicer! Now let us change the colours to black and red.

Grasping how programming languages handle colour can be a bit complex. The use of colour in data visualisation is an entire research field in itself. Randomly picking colours you think might go well together is a recipe for disaster. Lots of pre-made colour palettes exist for this reason. That said, we are just going to pick red and black! These are pretty contrasting and unlikely to cause confusion.

Define a colour palette in R

There are a number of ways of doing this, this is the simplest and embedded within ggplot without the need for another package.

cols <-  c("black","firebrick3", "purple", "green") # just a character vector, does nothing by itself, we only have 2 categories but we supply 4.

raw %>% ggplot(aes(x = Infection, y=`Myeloid cells`, col=Background)) + 
  geom_quasirandom(method = "tukeyDense",dodge.width = 0.75) +
  geom_boxplot() +
  scale_color_manual(values = cols) # function that looks up these colours

If we wanted to reuse the same palette in many plots we can extract the function and save it in R. Programming is all about being concise and anything that reduces copy and pasting will reduce error.

pal <-  scale_color_manual(values = cols)

raw %>% ggplot(aes(x = Infection, y=`Myeloid cells`, col=Background)) + 
  geom_quasirandom(method = "tukeyDense",dodge.width = 0.75) +
   geom_boxplot() +
  pal #  more consise

But what about the points being behind the boxplot - simple, reorder the layers!

raw %>% ggplot(aes(x = Infection, y=`Myeloid cells`, col=Background)) + 
  geom_boxplot() + # moved up a line
  geom_quasirandom(method = "tukeyDense",dodge.width = 0.75) +
  pal 

Editing thematic elements

We can change the axis titles, and just about anything else about the plot.

raw %>% ggplot(aes(x = Infection, y=`Myeloid cells`, col=Background)) + 
  geom_boxplot() + # moved up a line
  geom_quasirandom(method = "tukeyDense",dodge.width = 0.75) +
  xlab("") + 
  ggtitle("% Macrophages") +
  pal 

We can add bespoke themes such as cowplot using theme() function and its derivatives.

raw %>% ggplot(aes(x = Infection, y=`Myeloid cells`, col=Background)) + 
  geom_boxplot() + # moved up a line
  geom_quasirandom(method = "tukeyDense",dodge.width = 0.75) +
  xlab("") + 
  pal +
  theme_cowplot(font_size = 24)

Statistics

R is the preferred programming language for statistics because it comes inbuilt with lots of stats functions and was largely built by statisticians and ecologists who are basically statisticians.

We can run a t test for example using the t.test() function.

for_t <- raw %>%  filter(Infection == 'Infected')
t.test(`Myeloid cells` ~ Background, data = for_t)

    Welch Two Sample t-test

data:  Myeloid cells by Background
t = 12.235, df = 273.23, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 15.29692 21.16390
sample estimates:
mean in group C57BL/6  mean in group BALB/c 
             41.11024              22.87983 

Nice, but if we want to replace other graphing programs it would be nice to do this automatically on the plot. We can use the R package ggpubr to do this. Again we can store the function with specific arguments as an r object for future use.

comp <-   stat_compare_means(method = "wilcox.test", label = "p.signif", label.x.npc = "middle",size=8, vjust = 1,show.legend = FALSE)
t <-   theme_cowplot(font_size = 24)

raw %>% ggplot(aes(x = Infection, y=`Myeloid cells`, col=Background)) + 
  geom_boxplot() + # moved up a line
  geom_quasirandom(method = "tukeyDense",dodge.width = 0.75) +
  xlab("") + 
  pal + comp + t # much more concise 

Saving plots

Now that we have perfected our graph we need to save it. We save the plot as an object and print that to a file. We can do this with the ggsave function.

p <-raw %>% ggplot(aes(x = Infection, y=`Myeloid cells`, col=Background)) + 
  geom_boxplot() + # moved up a line
  geom_quasirandom(method = "tukeyDense",dodge.width = 0.75) +
  xlab("") + 
  pal +
  stat_compare_means(method = "wilcox.test", label = "p.signif", label.x.npc = "middle",size=8, vjust = 1,show.legend = FALSE) +
  theme_cowplot(font_size = 24)
ggsave(p, filename = "results/Myeloid_cells.pdf", width = 8, height = 5)

Advanced plots.

We have created a plot which is as high in quality as anything produced by prism or excel. R can create far more diverse plots than those programs, and you gain much more control over the visual aspects of the plots. However, it is a lot more complicated than using those GUI programs. The value of R (or Python) lies in not needing to manipulate the underlying data, e.g. by copy and pasting to create groups. Moreover we can perform transformations (functions) on the data simply, or change what we are looking out. We could even write a program that went through all your data and produced plots automatically - no more cumbersome copy and pasting.

We can perform maths on our data. Here we have been looking at % myeloid cells but maybe we want to look at total number of cells. \[\frac{ Cell Perentage}{100} * CellCount = Total number\]

We can do this in two ways:

1. Create a new col by multiplying Myeloid cells divided by 100 by by Cell Count using the mutate() function:

raw <- raw %>% mutate(`Myeloid cell number` = (`Myeloid cells`/100) * `Cell Count` )
raw %>% ggplot(aes(x = Infection, y=`Myeloid cell number`, col=Background)) + 
  geom_boxplot() + # moved up a line
  geom_quasirandom(method = "tukeyDense",dodge.width = 0.75) +
  xlab("") + 
  pal + comp + t

2. Perform the same maths within ggplot

raw %>% ggplot(aes(x = Infection, y=(`Myeloid cells`/100)*`Cell Count`, col=Background)) + 
  geom_boxplot() + # moved up a line
  geom_quasirandom(method = "tukeyDense",dodge.width = 0.75) +
  xlab("") + 
  pal + comp + t

Lets try to change the variable we are graphing. Remember you can print the variables using colnames(raw) function. ‘Live worms’ sounds interesting, let’s graph that, you can graph whatever you want.

raw %>% ggplot(aes(x = Infection, y=`Live worms`, col=Background)) + 
  geom_boxplot() + # moved up a line
  geom_quasirandom(method = "tukeyDense",dodge.width = 0.75) +
  xlab("") + 
  pal + comp + t

We can easily change the graph type to a scatterplot and add a linear regression

raw %>% ggplot(aes(x = Eosinophils, y=`Live worms`,col=Background, add = "reg.line")) + 
  geom_point() +
  geom_smooth(method = 'lm') +
  stat_cor(label.y = -2) +
  pal  + t + facet_wrap(~Background)

We can also look at change over time

raw <- raw %>% mutate(`Myeloid cell number` = `Myeloid cells` * `Cell Count` )
raw %>% ggplot(aes(x = Timepoint, y=`Myeloid cell number`, col=Background)) + 
  geom_point(position = position_jitterdodge(jitter.width = 10)) +
  geom_smooth() +
  pal + comp + t

ggsave(filename = "results/Myeloid_cells_timepoint.pdf", width = 8, height = 5)

Downlaod the data

So we will use a dataset from a paper I liked. This is a complicated paper but we want to take a take a small subset of the data and use it. This is for a macrophage in the peritoneal cavity called a large peritoneal macrophage (LPM). This is similar to the LCM in my dataset above.

We want to compare the transcriptome between LPM from WT mice given no stimulus and WT LPM given IL-4 complex a cytokine which activates type 2 immune responses and in macrophages is said to lead to ‘alternative activation’ (M2, M(IL-4) activation)

The paper has a link to the data in the methods section of the paper: GSE125727.

This brings us to a NCBI page for the datasets. For microarray datasets we can easily reanalyse these data using an R script built into the website called “Geo2R”. Analysis of RNA-sew data is a little more involved, but we will ignore that in this course.

I have highlighted the relevant samples. To start click geo2R.

On the next page, we define our groups. Do this by clicking ‘define groups’ typing 2 names, “naive” and “IL-4c”. To assign samples to those groups, cntrl or shift click the relevant sameps and then click on those groups.

How to assign groups

Then scroll down to the bottom of the page and click “Analyze”. This will tell the server to run the r script using the argument you provided (group information). You can actually view the R script if you like, or even copy it into R to run it on your own computer. Although you will need to install more packages, and the code might be a bit over your head at this stage. The analysis is a version of differential gene expression which looks at genes which are up or down in two groups.

This will take a couple of minutes to run, but eventually you will get a lot of outputs. What it will produce is a table of genes as rows, with various columns including 2 we are interested in, Fold change and P values.

Super cool outputs!

You can download these pictures if you like but what you really want is the table file. So select ‘download full table.’ This will download the a tsv file. A tsv file is like a csv file expect the data is separated by tabs not commas.

Move the file into the data folder for this project.

load in data

Here we will load in the data (assumming you were able to download it and put it in the right place). Unlike before we have to tell R that the tabs separate the data so the function needs another argument, sep="\t"

micro <- read.csv(file = "data/GSE125727.top.table.tsv", sep = "\t")

Lets investigate the file

dim(micro)

[1] 30177     8

So we have a lot or rows and 8 coloumns. Let’s look at the columns first.

colnames(micro)

[1] "ID"          "adj.P.Val"   "P.Value"     "t"           "B"          
[6] "logFC"       "Gene.symbol" "Gene.title" 

Hmmm, this looks fairly understandable. The cols we will be interested in are "logFC’ or log base 2 of fold change, Gene.symbol, and adjusted P value. We could select only that data to simplify matters.

micro <- micro %>% select(c("adj.P.Val","logFC","Gene.symbol"))

Now our dataset micro contains only three cols. The dataset contains 30177 rows, likely corresponding to genes (in reality they correspond to probes not genes, so the same gene may appear multiple times). We would not write out 30177 genes but we can use the head function to look at the first few.

head(micro$Gene.symbol)

[1] "AA467197" "Rnase2a"  "Chil4"    "Cst7"     "Pilra"    "Tmigd1"  

Cool! They look like genes to me! Chil4 is one of my favourite mouse genes too!

volcano plot

We will produce a basic volcano plot in R using ggplot2. We will need one more R package called ggrepel. Install this using install.packages("ggrepel").

library(ggrepel)

The most basic volcano plot can be made with the following code:

micro %>%
  ggplot(aes(x =`logFC`,
             y=-log(`adj.P.Val`) # do you notice the -log function?
               )) +
  geom_point() + theme_cowplot()

We can of course make this look better! Lets creat a new col that that can be used to highlight certain highly differential genes. The code from here on down might be a bit more complicated.

# mark up/down regulated genes
micro <- micro %>%
  mutate(gene_type = case_when(logFC >= 2.5 & adj.P.Val <= 0.05 ~ "up",
                               logFC <= -2.5 & adj.P.Val <= 0.05 ~ "down",
                               TRUE ~ "ns"))   

So what did that create? Let’s look

head(micro$gene_type)

[1] "down" "down" "down" "down" "up"   "down"

We will also create some variable that we can use for advanced plotting

# Add colour, size and alpha (transparency) to volcano plot --------------------
cols <- c("up" = "#ffad73", "down" = "#26b3ff", "ns" = "grey") 
sizes <- c("up" = 2, "down" = 2, "ns" = 1) 
alphas <- c("up" = 1, "down" = 1, "ns" = 0.4)

micro %>%
  ggplot(aes(x =`logFC`,
             y=-log(`adj.P.Val`),
             fill = gene_type,    
             size = gene_type,
             alpha = gene_type)) +
  ggtitle('naive LPM vs IL-4 complex LPM') +
  xlab('Log2 Fold Change') +
  geom_point(shape = 21, # Specify shape and colour as fixed local parameters    
             colour = "black") + 
  geom_hline(yintercept = -log10(0.001),
             linetype = "dashed") + 
  geom_vline(xintercept = c(-2, 2),
             linetype = "dashed") +
  scale_fill_manual(values = cols) + # Modify point colour
  scale_size_manual(values = sizes) + # Modify point size
  scale_alpha_manual(values = alphas) + # Modify point transparency
  scale_x_continuous(breaks = c(seq(-5, 5, 1)),  limits = c(-5, 5)) +
 theme_cowplot() +
  theme(legend.position = 'none')

That looks nicer. I won’t describe each line above, but you can read each line and see if you can work out what it does. You can comment out a line in the code using a hash # at the start of a line to stop the code being run, but without deleting it. Try commenting out various lines to see what that does.

Adding text to the volcano plot

The last thing we want to do is add text to the graph. Now we do not want to label all of the genes or our plot will look terrible and take hours to run. We just want to highlight the top few differential gene and highlight them on the graph.

First lets select the most differential genes. First arrange by LogFC ignoring minus sign using the abs() function and the desc() to place them in decending order.

micro <- micro %>% 
  arrange( desc(abs(logFC)))
head(micro, 20)

   adj.P.Val     logFC   Gene.symbol gene_type
1   9.67e-07 -8.391100       Rnase2a      down
2   2.28e-07 -7.661585      AA467197      down
3   1.07e-05 -6.095598          Mgl2      down
4   4.65e-06 -5.747173          Cst7      down
5   9.71e-06 -5.659000         Chil3      down
6   2.96e-06 -5.504030         Chil4      down
7   5.56e-06 -5.436653         Timd2      down
8   2.64e-05  5.359922          Fpr1        up
9   3.42e-05  5.320450        Cxcl13        up
10  4.95e-05 -5.071113          Ucp1      down
11  4.75e-05 -5.016402       Klk1b11      down
12  2.10e-05 -4.988765 I830127L07Rik      down
13  2.10e-05 -4.988538          Mall      down
14  3.81e-05 -4.849938        Gm5150      down
15  8.41e-06 -4.813952          Spp1      down
16  9.77e-06 -4.772435          Gatm      down
17  3.26e-05  4.520842         Cxcl2        up
18  9.10e-05  4.516250          Ccl3        up
 [ reached 'max' / getOption("max.print") -- omitted 2 rows ]

Now get the 32 top genes here.

genes_to_graph <- head(micro$Gene.symbol, 30)

What if we wanted to add some genes we like to this? Well we can just append the vector with some more gene symbols

genes_to_graph <- c(genes_to_graph, "Chil3", "Retnla", "Egr2", "Gapdh", "B2m" )

Next we want to create a new object which is a subset of micro that contains only the genes we are interested in. This is useful for telling ggplot to label only specific points. There are actually many ways of doing this and this is atually not a very intuitive method, but it works.

identify <- micro %>% filter(Gene.symbol %in% genes_to_graph  )

Finally we graph this.

micro %>%
  ggplot(aes(x =-`logFC`,
             y=-log(`adj.P.Val`),
             fill = gene_type,    
             size = gene_type,
             alpha = gene_type)) +
  ggtitle('naive LPM vs IL-4 complex LPM') +
  xlab('Log2 Fold Change') +
  geom_point(shape = 21, # Specify shape and colour as fixed local parameters    
             colour = "black") + 
  geom_hline(yintercept = -log10(0.01),
             linetype = "dashed") + 
  geom_vline(xintercept = c(-2, 2),
             linetype = "dashed") +
  scale_fill_manual(values = cols) + # Modify point colour
  scale_size_manual(values = sizes) + # Modify point size
  scale_alpha_manual(values = alphas) + # Modify point transparency
  scale_x_continuous(breaks = c(seq(-5, 5, 1)),  limits = c(-5, 5)) +
    geom_label_repel(data = identify,size=5, # Add labels last to appear as the top layer  
                   aes(label = Gene.symbol),
                   force = 3,
                   nudge_y = 1,max.overlaps = 200) +
 theme_cowplot(font_size = 18) +
  theme(legend.position = 'none')

ggsave(filename = "volcano_.png",width = 8,height = 5,dpi = 300)

Try to reverse engineer some of that code. I admit it is likely a bit too difficult for you.